#include<iostream>
using namespace std;
#include<cmath>

#include"library.h"

const int n = 100000;

double ratio( double R, int n ){
	
	int i, count = n;
	double r, costheta, phi, x, y, z;

	srand(time(NULL));
	for(i=0;i<n;i++){

		//①  取初始位置
		x = R; y = R; z = R;//这些值没有什么意义，只是方便进入下面的while循环
		while( x*x + y*y + z*z > R*R ){
			x = (double)rand()/RAND_MAX * 2 * R - R; // [-R,R]内均匀分布的随机数
			y = (double)rand()/RAND_MAX * 2 * R - R;
			z = (double)rand()/RAND_MAX * 2 * R - R;
			//取得的点在边长为2R的立方体内均匀分布
			//如果与原点距离大于R，会再次循环，重新取点
			//如果与原点距离小于R，则结束while循环，得到球体内均匀分布的随机点
		}
		//② 自由飞行
		r = (double)rand()/RAND_MAX * 2;//飞行距离为[0,2]内均匀分布的随机数
		// costheta在[-1,1]内均匀分布，phi在[0,2pi]均匀分布，这样取点，可以得到均匀分布的立体角
		costheta = (double)rand()/RAND_MAX * 2 -1;
		phi = (double)rand()/RAND_MAX * 2 * M_PI;
		x += r * sqrt( 1 - costheta * costheta ) * cos(phi);
		y += r * sqrt( 1 - costheta * costheta ) * sin(phi);
		z += r * costheta;
		// (x,y,z) + \vec{r} = 终点

		if( x*x + y*y + z*z <= R*R ){
		// 终点在球内，触发一次反应，得到2个自由中子，即自由中子数+1
			count ++;
		}
		else{
		// 终点在球外，消失不见，核弹内中子数-1
			count --;
		}
	}
	return (double)count/n;
}

double f(double R){
	return ratio(R, n) - 1;
}

int main(){

	double Rc = bisection(f, 0, 10, 1E-6);
	cout<<"The critical radius R = "<<Rc<<" (unit is mean free path of neutron)"<<endl;
	return 0;
}
